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Abstract 

The  present  research  effort  developed  a  real-space  formulation  for  orbital-free  density  func¬ 
tional  theory  and  Kohn-Sham  density  functional  theory  in  order  to  conduct  large-scale  elec¬ 
tronic  structure  calculations  that  take  an  important  step  towards  addressing  the  prevailing 
domain-size  and  geometry  limitations  of  existing  electronic-structure  codes.  In  particu¬ 
lar,  by  combining  the  real-space  formulation  with  a  finite-element  discretization,  which  has 
the  desirable  attributes  of  a  local  basis  that  is  amenable  to  coarse-graining  and  a  scalable 
discretization  on  parallel-computing  platforms,  it  has  been  demonstrated  that  large-scale 
electronic  structure  calculations  at  macroscopic  scales  become  accessible  to  conduct  an  ac¬ 
curate  electronic  structure  study  of  the  energetics  of  defects  in  materials.  Careful  verification 
and  validation  studies  of  the  developed  techniques  and  codes  for  both  orbital-free  DFT  and 
Kohn-Sham  DFT  have  been  conducted.  Using  orbital-free  DFT  the  energetics  of  vacancies 
in  Aluminum  as  well  as  the  properties  of  Al-Mg  alloys  have  been  investigated.  Using  Kohn- 
Sham  DFT  the  studies  on  the  energetics  of  Ni-Al  bi-layers  are  ongoing.  The  present  research 
study  was  conducted  in  close  collaboration  with  scientists  at  Army  Research  Labs.  In  this 
report  we  provide  a  summary  of  our  research  accomplishments  and  comment  on  our  ongoing 
and  future  work. 


1  Orbital-free  Density  Functional  Theory: 

Orbital-free  density  functional  theory  is  an  approximation  to  Kohn-Sham  density  func¬ 
tional  theory-which  is  widely  accepted  as  a  reliable  and  computationally  tractable  mate¬ 
rials  theory-where  the  kinetic  energy  of  non-interacting  electrons  is  explicitly  modeled  as 
a  functional  of  the  electron  density.  The  development  of  model  kinetic  energy  functionals 
for  orbital-free  DFT  is  an  active  area  of  research,  and  the  model  kinetic  energy  function¬ 
als  developed  over  the  past  decade  seek  to  capture  the  known  properties  of  uniform  elec¬ 
tron  gas  (Wang  &  Teter,  1992;  Wang  et  ah,  1998,  1999),  and  numerical  investigations  have 
demonstrated  the  accuracy  of  these  kinetic  energy  functionals  for  materials  systems  whose 
electronic  structure  is  close  to  a  free-electron  gas.  For  these  materials  systems,  orbital-free 
DFT  provides  a  computationally  efficient  approach  of  computing  material  properties  owing 


1 


to  the  linear-scaling  of  its  computational  complexity  with  systems  size  (as  opposed  to  a  cu¬ 
bic  scaling  for  conventional  Kohn-Sham  DFT  implementations).  The  main  aspects  of  the 
proposed  algorithmic  developments  in  orbital-free  DFT  include: 

(i)  Real-space  formulation:  The  ground-state  energy  of  a  materials  system  described 
by  DFT  is  given  by: 

E(p,  R)  =  Ts(p)  +  Exc(p )  +  Eh{p)  +  Eext{p ,  R)  +  EZZ{R )  (1) 

where  p  denotes  the  electron-density,  R  denotes  the  positions  of  nuclei,  Ts  is  the  ki¬ 
netic  energy  of  non-interacting  electrons,  Exc  denotes  the  exchange-correlation  energy, 

Eh  denotes  the  Hartee  energy  involving  the  Coulomb  interaction  between  electrons, 

Eext  denotes  the  electrostatic  interaction  between  electrons  and  ions,  and  Ezz  is  the 
nuclear- nuclear  repulsion  energy.  In  the  present  effort,  we  restricted  ourselves  to  LDA 
exchange-correlation  functionals,  and  developed  the  formulation  for  the  Thomas- Fermi- 
Weizsacker  family  of  orbital-free  kinetic  energy  functionals  as  well  as  the  Wang-Govind- 
Carter  (WGC)  kinetic  energy  functionals  which  involves  a  non-local  kernel  energy. 

The  various  terms  in  the  energy  functional  are  local,  excepting  the  electrostatic  in¬ 
teraction  and  the  kernel  energies  that  are  extended  in  the  real-space.  In  seeking  a 
local  real-space  formulation,  which  is  an  important  aspect  of  the  subsequently  devel¬ 
oped  coarse-graining  schemes,  we  developed  a  local  variational  reformulation  of  the 
extended  interaction.  In  particular,  the  extended  electrostatic  interactions  have  been 
reformulated  by  taking  recourse  to  the  solution  of  a  Poisson’s  equations.  Further, 
the  non-local  kernel  energies  have  been  reformulated  into  a  local  variational  formula¬ 
tion  by  solving  a  coupled  system  of  Helmholtz  equations.  The  details  and  the  associ¬ 
ated  mathematical  properties  of  our  local  reformulation  are  discussed  in  our  published 
manuscripts  Radhakrishnan  &  Gavini  (2010);  Motamarri  et  al.  (2012). 

(ii)  Finite-element  discretization:  A  general  finite-element  (FE)  framework  has  been 
developed  for  the  discretization  of  orbital-free  DFT  as  well  as  the  subsequently  dis¬ 
cussed  Kohn-Sham  DFT.  In  particular,  the  developed  framework  supports  higher-order 
finite-elements  up  to  8th  order,  including  spectral  finite-elements  which  are  better  con¬ 
ditioned  for  higher-order  discertizations,  ability  to  generate  coarse-grained  discretiza¬ 
tions  with  robust  meshing  scripts  for  CUBIT,  domain  decomposition  in  parallel  using 
PARMETIS.  As  part  of  this  research  effort,  we  investigated  the  computational  efficiency 
afforded  by  higher-order  finite-element  discretizations.  In  particular,  we  investigated 
the  convergence  rates  for  the  various  order  of  finite-element  approximations  and  the 
computational  efficiency  afforded  by  using  higher-order  discretization.  We  found  close 
to  optimal  rates  of  convergence  for  the  discretization  errors  despite  the  non-linear  na¬ 
ture  of  the  energy  functional  (cf.  Figure  1).  Further,  for  accuracies  commensurate  with 
chemical  accuracy,  we  found  significant  computational  advantage  by  using  higher-order 
basis  over  linear  FE  basis.  In  particular,  the  computational  advantage  was  over  100- 
fold  for  the  fourth-order  FE  basis  over  linear  FE  basis  (cf.  Figure  2),  with  diminishing 
returns  beyond  fourth-sixth  order. 

(iii)  Solution  strategies:  In  order  to  solve  the  non-linear  coupled  system  of  equations  de¬ 
scribing  the  electron-density,  electrostatic  potential  and  the  kernel  potentials  (c.f  Motamarri  et  al. 
(2012)  for  the  formulation),  two  strategies  have  been  explored:  (i)  A  staggered  solve 

where  for  any  change  in  electron-density  the  electrostatic  potential  are  recomputed; 
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Figure  1:  Convergence  rates  for  the  finite- 
element  approximation  of  bulk  Aluminum 
using  orbital-free  DFT  with  DD  kernel  en¬ 
ergy.  (Motamarri  et  al.,  2012) 


Figure  2:  Computational  efficiency  of  vari¬ 
ous  orders  of  finite-element  approximations. 
Case  study:  Aluminum  clusters  of  varying 
sizes.  (Motamarri  et  al.,  2012) 


(ii)  A  simultaneous  solve  where  all  the  coupled  systems  of  equations  are  solved  si¬ 
multaneously.  Our  mathematical  analysis  shows  that  the  staggered  solve  is  the  most 
robust  (Motamarri  et  ah,  2012),  though  it  is  computationally  also  expensive.  Our 
numerical  studies  suggest  that  the  robustness  of  the  simultaneous  solution  of  electron- 
density  and  electrostatic  potential  can  be  significantly  improved  by  using  improved 
preconditioners  like  block-Jacobi.  The  various  aspects  of  the  developed  solution  strate¬ 
gies  are  discussed  in  our  published  manuscript  Motamarri  et  al.  (2012). 

(iv)  Coarse-graining  using  quasi-continuum  reduction:  Building  on  the  ideas  pro¬ 
posed  in  Gavini  et  al.  (2007)  the  quasi-continuum  (QC)  reduction  for  orbital-free  DFT 
with  WGC  kinetic  energy  functionals  was  developed.  The  main  ideas  constituting  the 
development  of  this  coarse-graining  technique,  which  made  possible  orbital-free  DFT 
based  electronic  structure  calculations  on  multi-million  atom  system,  constitute:  (i)  the 
aforementioned  local  variational  reformulation  of  the  orbital-free  DFT  ground-state  en¬ 
ergy;  (ii)  an  effective  use  of  the  coarse-graining  ability  of  finite-element  basis  by  provid¬ 
ing  higher-resolution  where  necessary — e.g.  near  the  defect-cores — and  coarse-graining 
elsewhere.  Details  of  the  method  are  discussed  in  Gavini  et  al.  (2007)  for  TFW  kinetic 
energy  functionals  and  its  extension  to  account  for  WGC  kinetic  energy  functionals, 
developed  in  part  during  this  research  effort,  is  discussed  in  Radhakrishnan  &  Gavini 
(2010).  The  developed  QC  reduction  was  used  to  study  the  cell-size  effects  in  the 
formation  energies  of  vacancies  and  di- vacancies.  It  was  found  that,  even  for  simple 
defects  like  vacancies,  the  cell-size  effects  can  be  significant  up  to  1000  atoms.  Our 
studies  suggest  that  the  elastic  fields,  as  well  as,  the  perturbations  in  electronic  fields 
produced  by  defects  are  long-ranged,  and  that  cell-sizes  larger  than  those  used  conven¬ 
tionally  may  be  necessary  to  accurately  compute  the  energetics  of  defects.  Figure  3 
shows  the  cell-size  dependence  of  the  mono- vacancy  formation  energy. 

(v)  Transferability  studies  on  orbital-free  kinetic  energy  functionals  for  Al-Mg 
alloys:  As  part  of  the  present  effort,  a  study  was  also  conducted  to  validate  the  trans¬ 
ferability  of  the  orbital-free  WGC  kinetic  energy  functional  to  predicting  properties  of 
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Figure  3:  Cell-size  study  of  mono- vacancy  formation  energy.  (Radhakrishnan  &  Gavini,  2010) 


Al-Mg  alloys — a  material  system  of  significant  interest  in  developing  light-weight  ma¬ 
terials  for  transportation.  In  particular,  we  studied  the  phase  stability  of  these  alloys, 
and  the  results  are  presented  in  the  table  below.  The  /?'  alloy  has  879  sites  in  the  unit 
cell  with  a  disorder  in  20  of  these  sites.  In  our  studies  we  considered  the  two  extreme 
cases  of  all  the  20  sites  to  be  either  occupied  by  A1  (denoted  as  (3'(Al ))  or  by  Mg  (de¬ 
noted  as  P'(Mg)).  These  studies  show  that  the  orbital-free  kinetic  energy  functionals 
can  predict  the  relative  stability  of  the  various  phases  accurately,  which  is  a  stringent 
validation  test.  Further,  these  studies  also  demonstrate  the  efficiency  of  the  developed 
algorithms  which  allow  the  consideration  of  large  cell-sizes — for  instance  the  /?'  alloy 
containing  879  sites. 

Table  1:  Formation  energies  of  the  various  Al-Mg  alloys  computed  using  WGC  orbital-free 
kinetic  energy  functionals  using  bulk  local  pseudopotentials  (Zhou  et  ah,  2004).  Kohn-Sham 
DFT  calculations  using  ABINIT  with  Troullicr-Martins  pseudopotential  are  used  as  the 
validation  benchmark.  All  the  energies  reported  are  in  eV/atom,  and  negative  formation 
energy  denotes  an  energetically  favorable  alloy  formation.  Kohn-Sham  DFT  calculations  for 
the  f3 '  alloy  could  not  be  conducted  due  to  the  large  system  size. 


Alloy 

Al3Mg 

AluMg13 

Al^Mgn 

AI^qM  g23 

P\Al) 

P\Mg) 

Orbital-free  DFT 

-0.013 

0.069 

-0.007 

-0.001 

-0.026 

-0.084 

Kohn-Sham  DFT 

-0.009 

0.062 

-0.02 

-0.016 

NA 

NA 

2  Kohn-Sham  Density  Functional  Theory: 

While  orbital-free  DFT  is  attractive  due  to  the  linear-scaling  complexity  in  system  size,  the 
accuracy  of  the  kinetic  energy  functionals  limits  the  consideration  of  materials  systems  to 
those  whose  electronic  structure  is  close  to  free-electron  gas.  Thus,  as  part  of  this  research 
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Figure  4:  Computational  efficiency  of  vari¬ 
ous  orders  of  finite-element  approximations. 
Case  study:  Barium  2x2x2  BCC  cluster. 
(Motamarri  et  al.,  2013) 


Figure  5:  Computational  efficiency  of  var¬ 
ious  orders  of  finite-element  approxima¬ 
tions.  Case  study:  Methane  molecule. 
(Motamarri  et  al.,  2013) 


effort,  we  also  developed  computational  techniques  that  enable  large-scale  Kohn-Sham  DFT 
calculations.  The  ground-state  properties  in  Kohn-Sham  DFT  are  given  by  the  solution  of 
Kohn-Sham  equations: 


'H'lpi  =  CiV’i)  (2) 

where 

«=  (-iW  +  Vktp.R))  (3) 

is  a  Hcrmitian  operator  with  eigenvalues  eu  and  the  corresponding  orthonormal  eigenfunc¬ 
tions  for  i  =  1,2,  •••  ,N  denote  the  canonical  wavefunctions.  The  electron  density  in 
terms  of  the  canonical  wavefunctions  is  given  by 

N 

A>(r)  =X}l^(x)|2  ’  (4) 

i=  1 

and  the  effective  single-electron  potential  (Ves(p,R))  with  nuclear  positions  denoted  by  R 
in  (3)  is  given  by 

Ve&(p,  R )  =  Vext(R)  +  VH{p)  +  Vxc(p).  (5) 

In  the  above,  Vext,  Vh(p),  and  Vxc{p)  denote  the  potentials  associated  with  Eext,  Eh  and 
Exc.  Fourier  space  calculations  have  been  popular  for  Kohn-Sham  DFT  calculations  as  the 
computation  of  the  electrostatic  potentials,  which  are  extended  in  real-space,  can  be  conve¬ 
niently  computed  using  Fourier  transforms.  However,  a  Fourier-space  formulation  introduces 
periodicity  restrictions  which  severely  limits  the  study  of  defects  in  materials. 

In  the  present  work,  building  on  the  real-space  formulation  developed  for  orbital-free 
DFT,  we  implemented  a  real-space  formulation  for  the  Kohn-Sham  equations.  Furthermore, 
we  used  the  finite-element  discretization  of  the  formulation  which  enables  the  consideration 
of  complex  boundary  conditions,  also  provides  the  flexibility  to  use  an  adaptive  basis  set. 
Building  on  the  numerical  analysis  conducted  in  orbital-free  DFT,  we  computed  the  rates  of 
convergence  for  higher-order  finite-elements  for  both  pseudopotential  as  well  as  all-electron 
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Electron  Denslt 


10,03644 
0,02733 
0,01822 
0.009110 
0.0000 


Figure  6:  Electron  density  con¬ 
tours  of  7  x  7  x  7  FCC  aluminium 
cluster.  (Motamarri  et  al.,  2013) 


Electron  Density 


-31.96 

i-0.000 

Figure  7:  Electron  density  contours  of  100 
atoms  graphene  sheet.  (Motamarri  et  al., 
2013) 


calculations.  Further,  we  also  estimated  the  computational  efficiency  afforded  by  higher- 
order  finite-elements.  As  shown  in  figures  4  and  5,  for  chemical  accuracies  commensurate 
with  chemical  accuracy,  higher-order  finite-elements  provide  significant  computational  sav¬ 
ings  in  comparison  to  linear  finite-element,  discertizations.  Using  higher-order  finite-element 
discretizations  and  Chebyshev  acceleration  techniques  to  efficiently  compute  the  occupied 
eigenspace,  we  have  demonstrated  large-scale  electronic  structure  calculations  for  both  pseu¬ 
dopotentials  and  all-electrons — for  instance  a  7  x  7  x  7  FCC  cluster  containing  1688  atoms 
(pseudopotential  calculation)  and  a  graphene  sheet  containing  100  atoms  (all-electron). 

Computations  on  reactive  nanofilms:  Building  on  our  developments  in  Kohn-Sham 
DFT,  we  incorporated  non-local  Troullier-Martins  pseudopotential  in  the  Kleinman  Bylander 
form  and  conducted  studies  on  Ni-Al  bi- layers.  In  particular,  we  computed  the  excess  energy 
in  the  bi-layer  system  as  a  function  of  the  bi-layer  thickness.  This  excess  energy  is  computed 
as  the  difference  in  the  energy  of  the  Ni-Al  bilayer  system  and  the  energy  of  bulk  Al  and 
Ni  atoms  with  similar  composition.  Figure  8  shows  the  excess  energies  of  a  Ni-Al  bilayers 
along  (100)  as  a  function  of  the  bi-layer  thickness.  For  small  bi-layer  thickness  this  excess 
energy  is  negative  suggesting  that  the  decrease  in  the  energy  due  to  bonding  at  the  interface 
is  greater  than  the  increase  in  the  energy  due  to  lattice  mismatch.  However,  with  increasing 
bilayer  thickness  this  excess  energy  increases  proportionally  with  the  bi-layer  thickness,  which 
is  a  result  of  the  elastic  energy  from  the  lattice  mismatch.  Beyond  a  particular  bilayer 
thickness  it  becomes  energetically  favorable  to  nucleate  threading  dislocations.  In  our  future 
investigations  we  seek  to  compute  the  critical  bilayer  thickness  at  which  the  dislocations 
nucleate. 


3  Development  of  Software: 

In  collaboration  with  a  team  of  scientists  at  Army  Research  Labs,  which  includes  Mr.  Ken¬ 
neth  Leiter,  Mr.  Joshua  Crone,  Dr.  Michael  Scott,  Dr.  Jaroslaw  Knap  and  Dr.  Peter 
Chung,  we  have  committed  significant  time  and  effort  towards  developing  a  software  out 
of  the  computational  techniques  developed  over  the  course  of  this  project.  This  software 
is  being  jointly  developed  by  the  group  members  of  the  PI  and  the  team  from  ARL,  and 
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Figure  8: 


Excess  energy  as  a  function  of  bi-layer  thickness  of  Ni-Al  bi-layers. 


is  maintained  in  a  repository.  There  is  an  active  ongoing  collaboration  with  ARL  in  the 
development  of  this  software  with  weekly  conference  calls  and  annual  visits  by  PI  and  his 
students  to  ARL. 
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(revision  submitted,  2012). 
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(2012). 
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